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In-plane dissipation as a possible synchronization mechanism for terahertz radiation from intrinsic 

Josephson junctions of layered superconductors 
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Strong terahertz radiation from mesa structure of P^S^CaC^Og+a single crystal has been observed recently, 
where the mesa intrinsically forms a cavity. For a thick mesa of large number of junctions, there are many 
cavity modes with different wave vectors along the c-axis corresponding to almost the same bias voltages. The 
mechanism responsible for exciting the uniform mode which radiates coherent terahertz waves in experiments is 
unknown. In this work, we show that the in-plane dissipation selects the uniform mode. For perturbations with 
non-zero wave numbers along the c-axis, the in-plane dissipations are significantly enhanced, which prevent the 
excitation of corresponding cavity modes. Our analytical results are confirmed by numerical simulations. 

PACS numbers: 74.50.+r, 74.25.Gz, 85.25.Cp 



I. INTRODUCTION 

A layered cuprate superconductor, such as 
Bi2Sr2CaCu20s + 5 (BSCCO), intrinsically forms a stack 
of Josephson junctions 1 . Because of the large superconduct- 
ing energy gap (60 meV), these build-in intrinsic Josephson 
junctions (IJJs) can support oscillations with frequencies in 
the terahertz (THz) band. IJJs are homogeneous and packed 
on nanometer scale, much smaller than THz electromagnetic 
(EM) wavelength. If synchronized, they can radiate powerful 
THz EM waves. The radiation frequency is determined by 
the bias voltage according to the ac Josephson relation, thus 
in principle can be tuned continuously. The THz generator 
based on IJJs thus is promising to fill the THz gap 2 3 . 

In 2007, coherent radiations from a mesa structure of 
BSCCO without external magnetic fields were detected 
experimentally 4 . The radiation frequency and voltage fol- 
low the ac Josephson relation, thus the radiation is due to the 
Josephson plasma oscillation in the mesa. The frequency / 
is determined by lateral size L x of mesa / = co/(2L x ) with 
c o - cj y[e~ c the Josephson plasma velocity where e c is the di- 
electric constant of BSCCO. The measured relation between 
L x and / has revealed unambiguously that the mesa works as 
a cavity to synchronize plasma oscillations in different junc- 
tions. The cavity resonance mechanism has been confirmed 
by many independent experiments^. 

The experiments raise several questions. First in experi- 
ments a dc current is uniformly injected into the mesa. One 
thus expects the superconducting phase would oscillate ho- 
mogeneously along the lateral directions, which, however, 
seemed to be hardly reconciled with the observed cavity 
modes. This question has been addressed in Refs. [TOland 
ITTI They suggested that superconducting phase develops n 
phase kinks near the cavity resonances. With the help of the 
jt phase kink, the standing wave of electromagnetic fields can 
be stabilized and a large amount of energy is pumped into EM 
waves from the dc current. 

Secondly, as the mesa intrinsically forms a three- 
dimensional (3D) cavity, cavity modes both along the c- 



axis and ab-plane can be excited. When the frequency 
of the plasma oscillation co determined by the bias volt- 
age per junction is tuned to the cavity frequency a>' c = 
■\J(m x 7T/L x ) 2 + (m y jr/Ly) 2 Cq, the cavity mode (m x ,m y ) is 
excited 12 , where c q depends on the wave vector q - nn/(N+l) 
along the c-axis with N the number of junctions and n an 
integer. For a large N, such as N ~ 1000 in experiments, 
there are many cavity modes with different q within a narrow 
window of bias voltage, and one would expect various modes 
should be excited when the voltage is swept. Nevertheless, for 
a given radiating sample, only the modes uniform along the c- 
axis (q - 0) is observed upon sweeping current in experiment. 
The reason remains illusive. 



Thirdly, understanding on possible cavity modes along the 
c-axis is also important for getting stronger radiation from the 
mesa structure of BSCCO, which is still too weak for practical 
applications to date. The radiation power can be enhanced by 
using thicker mesa with larger N, since the radiation power is 
proportional to N 2 in the superradiation region. One question 
is whether there exists a fundamental limiting factor besides 
the heating effect for N. 



In this paper, we show analytically that in-plane dissi- 
pations prevent the excitation of non-uniform cavity modes 
along the c-axis. For non-uniform perturbations with a large q, 
effective in-plane dissipations are greatly enhanced, and thus 
the perturbations quickly die out and no cavity mode is ex- 
cited. For perturbations with a small q, in-plane dissipations 
are weak and vanish for q - 0. The weak dissipation along 
the c-axis cannot damp them efficiently, and thus cavity modes 
with small g's are excited. For N « 10 3 , only the cavity mode 
with q - can be excited. While for large N « 10 4 , modes 
with finite but small q vectors can also be excited and com- 
pete with mode q - 0. The analytical results are confirmed by 
direct numerical calculations. The mechanism to achieve syn- 
chronization by dissipation may be applied to other systems 
as well by preventing the excitation of the out-of -phase mode. 



II. MODEL 

We consider a stack of IJJs with lateral sizes L x a 80 /im, 
L, ~ 300 yum similar to those in experiments, see the inset of 
Fig. fl] Because L y » L x , we can assume that the supercon- 
ducting phase is uniform along the y-axis. The dynamics of 
the gauge invariant phase difference <pi and magnetic field B Y j 
in the /-th junction are described bySEl 



d t <pi +/3 c d,(pi + sinipi = d x B yJ , 



(1) 



[<TA (2) - ( 1 + B ab d,)\ B y j + ( 1 + B, b d,) djpi = 0, (2) 

where A (2) // = fi+i + //_i - 2fi is the finite difference op- 
erator. ( x 10 5 is the inductive coupling, B c ss 0.02 and 
B a i, =» 0.2 are the renormalized conductivity along the c 
axis and ab plane respectively. 18 For BSCCO, it is well es- 
tablished that the in-plane conductivity cr a t, is much larger 
than the c-axis conductivity cr c , from microwave 19 , infrared 
spectroscopy 20 and transport 21 measurements. At tempera- 
ture T = 0, (T ab ~ 4 x 10 6 (Q • m)- 1 and lt c ~ 0.2 (Q • m) _1 . 
Both B c and B a \, depend on frequency, but t he dependence is 
weak in the interested frequency region 22 23 . In the following 
discussion we will neglect the frequency dependence. The 
generalization is straightforward as we are working in the fre- 
quency domain. The in-plane dissipation due to B a b has been 
overlooked in many theoretical models. 

Equations (fill and Q are supplemented by the boundary 
conditions. When the phase oscillates uniformly along the 
c-axis, strong radiation of EM waves occurs, which can be 
accounted for using the boundary conditio n 1 24 ! 25 ! 
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Z(o>r ' ys^j-^hilf]' 



(3) 



where - (+) corresponds to the edge x = L x (x - 0), and k^ - 
to yfej with €d the dielectric constant of the dielectric medium 
outside the IJJs. For stacks with height L z <k 100 p.m, Z » 1. 
For non-uniform oscillations along the c-axis, the radiation is 
weak and we can use the non-radiating boundary condition 
B y j = +I ext L x /2 with 7 ext the bias current. We assume that the 
IJJs stack is sandwiched by two good conductors, such that 
the tangential current inside the conductor is zero. We use the 
boundary condition B v ./=i - B v j = o and similarly for I = N, 
which corresponds to d z B y (z) — in the continuum limit. 



III. INSTABILITY OF THE HOMOGENEOUS SOLUTION 

In experiments, one first ramps up the bias current and when 
current exceeds the critical one, the IJJs switch into the resis- 
tive state. One then reduces the current to the target value 
where radiation is observed. In the resistive state with the 
current close to the critical one, the phase tpi oscillates homo- 
geneous along the x direction. The phases may be either uni- 
form along the c-axis or different in different junctions. Here 
we show that the solution with <p\ homogeneous along the x 
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FIG. 1. (color online) Stability diagram for the homogeneous solu- 
tion along the ab plane. In the filled region, the homogeneous so- 
lution is unstable and the cavity mode ik m ,q) is excited. Inset is a 
schematic view of the setup for THz radiation. 



direction is unstable when the bias current is reduced and that 
only the cavity modes with a long wavelength ^« 1 along 
the c-axis can be excited. 

For the uniform solution along the c-axis tpi = tpo, the solu- 
tion to Eqs. (fTJ and (J2]i in the THz frequency region a> » 1 
can be written as cpo - a>t + g(x) exp(z'wf) with 



g(x) = 



1 



+ iB c <jj 



cos [(x - L x /2)oj] 



(4) 



Z sin (L x a>/2) - i cos (L x a>/2) 

The first term in the square bracket of Eq. Q is due to plasma 
oscillation and the second term is due to radiation. The fre- 
quency co can be tuned by the bias current 7 ext (or voltage). 

To reveal the instability of the homogeneous solution, we 
add small perturbations to the uniform solution tpi - <fo + 
6i, B y i - B + B yl with B - d x ip , 9 ( «1 and B y j <K 1. 
Since the perturbations are non-uniform along the c-axis, the 
radiation contribution can be neglected, thus we can use the 
non-radiating boundary condition, d x 0i - and B y j = 0. The 
solution for 6/ and B y j can be written as 

8i(x, t) — y cos(ql) cos (k m x) a p (m, q) cxp[i(pco — D.)t] (5) 



B y j(x,t)- y cos(ql) sin (k m x) b p (m, q)e\p[i(pco — D)t] (6) 

m,q,p 

with k m - mn/L x , q - nn/(N +1) and p an integer. The per- 
turbations with frequency Q. couple with the nonlinear Joseph- 
son current sin if/ and induce frequency harmonics pco — Qpy 
From Eq. Oh, we have b p - -c 2 q k m a p with the plasma veloc- 
ity c q (co) - [1 +2£(1 -cos q)/(l +B ab ico)Y xil . Substituting 
Eqs. p) and (J6l) into Eqs. (fTh and d2), and comparing each 
frequency component, we obtain equation for perturbations 



[t<4(Wp) - co 2 p + icOpBc] a p 



cip-i + flp+i A 



p-2 



2/ 



= 



where u? m - c 2 fc 2 , u) p - pa> - Q,, and A p - 
Y L " dxg(x)9p(x) cos(k m x). Because of the in-plane damp- 
ing p a b, w 2 is a complex number. We can split a> 2 n into real 
and imaginary parts oJ^co) = Q 2 (a>) + i!B ab (a))oj, with 



"„ 2 » 



l+2^(l-cos g )+^ 2 . 
[1 + 2f(l- cos q)] 2 +pl h aj 2 ' 



5ab 



(«) 



2^(1 -cosq){3. db 
[l+2£(l-cos (? )] 2 + > Si 



> 2 



(7) 



(8) 



£l m (co) is the cavity resonance frequency for the mode (k m , q). 
For the non-uniform perturbations along the c-axis q > 0, the 
effective in-plane damping is enhanced according to < B db . It is 
this enhanced in-plane dissipation that prevents the excitation 
of non-uniform cavity mode along the c-axis with a large q as 
revealed later. 

In the region of u) m » 1 and j3 c «: 1, we have |Q 2 -w 2 () | <K 
1. Thus the dominant wave vector of p (x) is k m . The dom- 
inant wave vector for g(x) is k x - because the radiation 
contribution is small. Then we can approximate 8 p (x)g(x) by 
g6 p (k m ) cos(k„,x) with g the spatial average of g(x), because 
other modes are negligibly small. Neglecting the small dissi- 
pation contribution /3 C in g, we have 



1 



1 



1 



2i -2w 2 L x u) 3 [cot (L x u>/2) + Zi] -2oo 2 



R r + iRj 



where R r is the real part of the radiation contribution and R, is 
the imaginary part. We then have A p - ga p . The equation for 
perturbations then can be written as 



^m - w ;; + R r ~ 1T~ j 



a p + i \o) p (j3 c + S ab ) + Ri] a F 



1 / x 1 

+ 2 \ a i'- 1 + a P +x ) + Y~i ap ~ 2 ~ RrG i'- 2 ~ iR i a P-2 ~ °- ( 9 ) 

The radiation shifts the resonance frequency by R r and also 
contribution to the damping through /?,-. Both R r and Ri have 
the order 1 /a> 2 <K 1 for L x ~ 1 . The resonant frequency for 
perturbations with a wave vector (k m ,q) is Q. — a> c with a> c 
given by a> 2 = Q 2 (<y c ) + R r ((L>) — ^, u> c is also the cavity 
frequency for the cavity mode (k„,,q). 

We then show that the parametric instability develops at a 
voltage when to - 2u c + 5 with 5 <K 1 . As Q » Q. m w co c , the 
dominant frequency components are p — and 1. The other 
frequency harmonics (p > 1 or p < 0) are small because their 

amplitude is of order ([(2/7 - l) 2 - l]Q 2 ,j flo.i <s «o,i thus 
can be neglected. Then Eq. (|9|l can be written as 

[-2io c £l s + i (-o) c (J3 C + S. db ) + Rd] «o + ai/2 = 0, (10) 



The homogeneous solution along the x direction is unstable 
when Im[fi,5] > 0, which gives 



5 2 < 



(fid,(o» c )+A) 



1 

2^? 



1 



rff 



w- 



(12) 



From this expression, it becomes clear that the radiation tends 
to destabilize the homogeneous solution, albeit with a small 
impact Ri ~ 1 /u> 2 <K 1 . More importantly, both the effective 
in-plane dissipation S a b and dissipation along the c-axis /3 C 
tend to stabilize the solution homogeneous along the x axis. 
For long wavelength perturbations along the c-axis, £q 2 <k 1, 
the dissipation is weak S a b (o> c ) +/3 C < l/(2w c ) thus the cavity 
mode (k m , q) can be excited. By using Eqs. d7), dS} and ( 12 1, 
the condition for the excitation of the mode (k m , q) for /3 a t, > 
is 

^ < 4( A "i)( 1+ ^^ +A "i) 1 - (13) 

The stability diagram is presented in Fig. [T] with the given 
parameters. Only the cavity mode with long wavelength along 
the c-axis can be excited. As q - rm/(N+l), for N ~ 10 3 used 
in experiments, only the uniform cavity n — can be excited. 
For a larger N > 10 4 , modes with finite but small n can also 
be excited. 

For the solution that is homogeneous along the jt-axis but 
non-uniform along the c-axis, the radiation contribution can 
be neglected and the phase oscillates according to ipi = <f>i + 
cat + i/{-u> 2 + ijioo), where (pi accounts for the phase shifts in 
different junctions and is randomly distributed. The stability 
analysis is the same as the case of uniform solution, but now 
there is no radiation contribution R r = R, - 0. The stability 
diagram is the same as that in Fig. [T] because the radiation 
contribution is small for L x « 80 fim used in experiments. 

For a small crystal proposed in Ref. [25] L x » 4/vm, 
oj c ~ 150. The cavity modes cannot be excited for such a 



high frequency according to Eq. ( 13 I. The homogeneous so 



lution along the lateral directions is thus stable. In the single 
junction limit N — 1, the parametric instability leads to the ex- 
citation of solitons, which manifests as the zero-field current 
step in IV characteristics^. 

For non-uniform perturbations along the c axis, the in- 
plane current J x j is induced according to the Ampere's law 
4nJ x j/c - -(B,.,/+i -Byj)/s where we have neglected the 
displacement currenfSI. J x j has contribution from the normal 



current and supercurrent J x j - cr^jr.dt'Pi + „ 2 % Pi where 



c4>0 



8^:, 



Pi = d x &i - 2nA x /<f>o is the in-plane superconducting mo- 
mentum with superconducting phase 0/ and vector potential 

A x . The in-plane dissipation in units of ^ (^Tks) ^ s 



[-2w c (6 - Qj) + i (oo c Q3 C + S, b ) + /?,)] a x + a /2 = 0, (11) 

with Q = Q.s + u> c and Q^ «c 1 . The spectrum of the perturba- 
tions CI is given by equating the determinant of the coefficients 



matrix in Eqs. (|T0|> and ( 1 1 1 to zero, which yields 



Cls 



2 



2 2 



l l r, I 

— - -oj c 6 2 - 6Ri + —R] 

4oj r 2 2(jJ r 



Pab(km, q, (O) = 



2sm 2 (q/2)/3. db OJ 2 
(AbOj) 2 + 1 



\4\^ m \a p (k m ,qf , (14) 



which is much larger than the dissipation along the 
c-axis P c (k m ,q,oo) - \oj 2 fi c \a p (k m ,q)\ for g's when 
sin 2 ( ? /2)>[0S ab w) 2 + l]/[4^ b |c 4 |4]. Thus the non- 
uniform perturbations with a large q quickly die out due to 



(a) I„ t = 0.27, ft, = (b)/ cxt =0.30, ft» = (c)/„ t = 0.57 

rif *ff ifit 




= (d) l ex 



3.49 ft, = 



IIIIHIMikiiiin|' V ^ 



£ 200 



IMMMIH 



ifiitti 

••••••• 

••••••• 



u 



Iffllvflfl 

IIMlllll 
MIIUMI 

Hi. A. lit 







le) I ezt = 0.16, ft 6 = 0.2, M = 10 3 (f) I ea = 0.18, ft„ = 0.2, W = 10" 




0.2 0.4 

X[A C 



0.4 0.2 




FIG. 2. (color online) Snapshots of the electric field (first row), magnetic field (second row) and Josephson current sin(</)/) (third row) without 
the in-plane dissipation/?,^ = (a-d) and with the in-plane dissipation fi ab = 0.2 (e-f). When/3,,;, = 0, various cavity modes (m,n) - (7,5) in 
(a), (m,n) = (9,3) in (b), (m, ri) = (7,2) in (c) and (m, ri) = (4, 1) in (d) are excited when the bias current is swept. For/? n6 = 0.2, only the 
modes with n = 0, [(m, ri) = (1, 0) in (e), other cavity modes with m > 1 are not shown here] are excited for N = 10 3 . For N = 10 4 , irregular 
patterns of the EM fields are excited even with the in-plane dissipation/?,^ = 0.2. The supercurrent forms blocks with alternating sign between 
neighbouring blocks indicating ±n phase jumps at the interface of each block. We use L x = 0.4A C , f = 7.1 x 10 4 in simulations. 



the strong in-plane dissipation P a \,. While for perturbations 
with g«l, the in-plane dissipation is weak or absent, and the 
perturbations lead to the excitation of the cavity mode with a 
small q. 



IV. NUMERICAL SIMULATION 

We also perform numerical simulations by solving Eqs. 
(fill and dZk with the non-radiation boundary condition B y j - 




20 30 40 

Voltage co [<1> co /(2rcc)] 

FIG. 3. (color online) IV curve obtained with/? ai = 0.2 and N = 10 3 . 
Other parameters are the same as those in Fig. [2] The radiation 
power is estimated using 5 r = |£„,.|/(2|Z|) and e rf = 1 where E ac is 
the averaged ac electric field at edges 28 . The maximal power at the 
first cavity mode in = 1 is about 2800W/cm 2 . The experimental 
measured IV deviates significantly from the theoretical one in high 
bias region due to the strong self-heating effect. 



±I ext L x /2 to check the above analytical results. For numerical 
details, see Appendix A. Let us first present the results with- 
out in-plane dissipation/?^ = as shown in Fig. |2Ta-d). Upon 
sweeping the current, various cavity modes both along the x- 
axis and c-axis are excited for Af = 10 3 . It is difficult to excite 
the uniform mode q — due to the existence of other com- 
peting cavity modes in the crystal with this large Af at a given 
voltage. Interestingly, when we turn on the in-plane dissipa- 
tion by putting /? fl fe = 0.2, only the cavity modes uniform along 
the c-axis can be excited when the current is swept as shown 
in Fig. |2|e), consistent with the above analytical results. The 
corresponding IV characteristics is shown in Fig. [3] At the 
cavity resonances, a large amount of energy is pumped into 
the plasma oscillation and current steps are induced, and the 
radiation power is enhanced. 

For N < 4000 in simulations, only the uniform mode can 
be excited due to the in-plane dissipation, consistent with the 
results in Fig. [T] For even larger number of IJJs, such as 
N - 10 4 , irregular patterns of EM fields are developed in sim- 
ulations for the adopted parameters, see Fig. f2\f). It is consis- 
tent with the analytic result summarized in Fig. [Tlwhere cav- 
ity modes with n - 0, 1, 2 can be excited. Other mechanisms 
to synchronize all junctions are needed in order to achieve 
strong radiation for thick mesas. Here, we note that for small 
values of N < 100, the voltage corresponding to various cav- 
ity modes is discrete because c q is well separated for different 
values of «, thus one can select the cavity modes simply by 
tuning the voltage. 



For p a i, - 1.0, we found numerically that the uniform 
plasma oscillation with q — becomes stable again for N = 
10 4 . This suggests that stronger in-plane dissipation is impor- 
tant to achieve uniform plasma oscillation. B a b increases with 
temperature, while on the other hand thermal fluctuations tend 
to destroy the uniform oscillation. Therefore the synchroniza- 
tion is optimal at some intermediate temperatures. 

One peculiar feature in Fig. [2] is that supercurrent forms 
blocks in space, where the current change sign between neigh- 
boring blocks, or equivalently at nodes of the oscillating elec- 
tric field. This means that there is +n phase jump or ±n phase 
kink at the interface of blocks. The phase kinks stack along the 
c-axis with alternating signs, such as (■ • ■ , 1, — 1, 1, — 1, ■ • ■ )tt. 
The state with phase kink with q — was first suggested in 
Refs. [TOlandlTTIas a possible mechanism for strong THz radi- 
ation observed in experiments. A characterization of the kink 
sate with a general q is presented in Appendix B. 



V. DISCUSSIONS 

One should also check the stability of the excited cavity 
modes. This has already been done in Refs. 28-30 and found 
that the kink states associated with cavity modes with q = 
are stable for a small N, which is also confirmed by the 
results in Fig. [2] For a large N, it was shown in Ref. 30 that 
long wavelength instability develops and the kink states with 
uniform plasma oscillation q - becomes unstable. This is 
consistent with the results in Fig. [2jf). 

When a strong magnetic field is applied parallel to the ab- 
plane, the Josephson vortices (JVs) are induced. The JVs 
favor the triangular lattice in low velocity region due to the 
strong inter- vortex repulsion in a long IJJs stack 31 . As shown 
in Appendix C, our simulations show that the in-plane dissi- 
pation mechanism becomes insufficient to achieve the rectan- 
gular JVs lattice and in-phase plasma oscillation, consistent 
with the analytical results in Ref. [TTl A new mechanism for 
synchronization is needed in this case. 
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Appendix A: Numerical method 

Here we present the numerical method to solve Eqs. (fill 
and pi. The model is discrete inherently along the c-axis, 
and we only need to discretize along the x-axis. The phase 
<p is defined at nodes (J, I) and magnetic field B y is defined at 
nodes (j + 1/2, /) with integer j and I, see Fig. |4] The grid size 



along x is dx and time step is dt. Equation (fT) then becomes 

*>• ' ■ m 

+ sin ipjj 
B" 



dfi Pc 



2dt 
E 



IJ+1/2 



'U-l/2 



dx 



(Al) 



We know <p and B y at the m-th time step and we can obtain 



<p at the (m + l)-th step directly from Eq. ( Al I. We use an 



implicit method to discretize Eq. |2|. After some simple ma- 
nipulations, we have 
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<PlJ-\ + <Pij+i - fij-i 



2dx 
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. «i-<i)-kfri-y&-i) 

ab T~, 



dtdx 
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XM/2 

dt 



(WA«) 



'lj+1/2 



nm+l Dm+1 

a ( WA «)^ +Ab W 



(A2) 



Equation ( |A2| > can be written as a matrix equ ation. Inverting 
the matrix at the right-hand side of Eq. (A2i using ip'] l+l ob- 



tained from Eq. (|A1||, we then obtain B'" + ^ l/2 . The electric 



field is given by Ef. = « +1 - ip" l 7 l )/(2dt). 

h] HJ hj 

Appendix B: Characterization of the kink state with a general q 

In the kink state, (pi can be written as 
(pi - cot + <p s j(x) - i 2_j A m, q cos(k m x) cos(ql) exp(icot), (Bl) 

m.cj 

where the rotating phase at the right-hand side (rhs) of Eq. 



( B 1 1 is due to the voltage, the second term at rhs is the static 
phase kink and the last term at rhs is the cavity mode both 
along the c-axis and x-axis. For stacks with a large N, there 
are many cavity modes {m,q) corresponding to a same fre- 
quency co, thus a summation over all modes are needed 28 . The 
magnetic field is 

Bi - B s j(x) - i } C mA sin(£ m x) cos(gf) exp(icot), (B2) 

m.q 

where B s j is the static magnetic field due to the phase kink. 



We have neglected the frequency harmonics in Eqs. (Bl i and 



z 1 + 2 
l + l 

I 
1-1 



V&2J-1 






T 



1+2J-1/2 



<Pl + 2J 



T 



r 



1+2J+1/2 



Vu2,i+\ 



Piiu-i 






"i+lj-1/2 [: fi+ij \ D i+i,j+i/2 



-1/2 



Vu 



Bl 



J+l/2 



"Pl+lj'+l 

~l 



•Pij+i 



i-i J- 1/2 

>x 



J 

i + 1/2 i + 1 



FIG. 4. (color online) Schematic view of the numerical grids. 
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FIG. 5. (color online) Snapshots of the Josephson current sin(<^/) in the flux-flow region with the in-plane dissipation p ab = 0.2. Here the 
applied magnetic field is B a = 1 T, L x = 0.1/t c , N = 500, f = 7.1 x 10 4 . The core of Josephson vortex is located at <p = (2m + l)n. 



c, 



B2| >, which is valid when A„ hq < 1. From Eq. (J2l we obtain 
A m ,qk ln c 2 and [1 - £A {2) ]B sJ (x) = d x tp s j(x). Substi- 



/n.ty 



tuting Eqs. (Bl i and ( |B2| i into Eq. (fill, we obtain a closed 
equation for </?/. For the frequency component a>, we have 

lik^c 2 - (b c (jj + ico 2 j\ A„ uq cos (fc m jt) cos (ql) - -ie mj . (B3) 

Projecting exp(/<^ s ) into the cavity mode (m, q), we obtain the 
amplitude of the plasma oscillation 



— iF 

. ll m,q 

m ' q iklf 1 , - (J3 c u> + ico 2 ) ' 



(B4) 



with the coupling between the cavity mode and phase kink 



a N r 



exp(iip s j) cos (k„,x) cos (ql) dx, (B5) 



where a - 2 for q — and a - A for q > 0. When the voltage 
is tuned close to the cavity resonance, the amplitude of plasma 
oscillation is enhanced. The linewidth of the resonance is de- 
termined by p a b and f3 c . The TV characteristics is given by 
4xt = fic<J-> + ( sm (<^/))x,/,/ where (• • • ) xlt denotes average over 
space and time. We then obtain 



4xt = y6 c W + Re 



2-< ik 



\F m , q \ 2 /(2a) 



$L& - (B c (o + ioj 2 ) 

m,q m q \r-t / 



(B6) 



Now let us consider the static component 
%p M (x,z) = yA (2) Yj [A m , ? cos(fc m x)cos( ? /)e-^'] . 



(B7) 



The solution of ip s j depends on the spatial profile of the 
plasma oscillation. To present the analytical results, we con- 
sider a case that the plasma profile has well-defined nodes as 
shown in Fig. I3](a-d), where we may approximate the spatial 
profile by a dominant mode Acos(&ix)cos(g/). Without loss 
of generality, we have taken the m = 1 mode. The variation 
of ip s j along the c-axis is <p s j - (— l)V*o except for the node 
region of the plasma oscillation, which is much faster than the 
plasma mode q. We may approximate Eq. (|B7|> as 



d x (pso = 2£Re[A]cos(&ix)cos(gf) sm ( l ^o)- (B8) 



Equation ( |B8| > is invariant under the transformation x <— 
L x — x and tp s o <— n - ip s o, which gives the static n phase 
kink along the jc-axis at the nodes of oscillating electric 
field when cos(kix) = 0. The width of the kink is A^ = 
1 / J2£Rs\A]\ cos(ql)\, which is small 4« 1 except for the 
node region of cos(ql) » 0. The kink in the Z-th junction can 
be excited only when A^ <K L x . Near the nodes of the oscil- 
lating E z or By along the c-axis where cos(ql) ss 0, Ak may 
be comparable to L x thus no kink exists in the node region, 
consistent with results in Fig. [2] When cos(ql) changes sign, 
(fi S Q acquires a n shift, because Eq. (B8i is invariant when 
cos(ql) <— -cos(ql) and ip s o <— n + ip S Q. Thus there are ±n 
phase jumps at the nodes of oscillating electric fields both 
along the c and x axis. 



Appendix C: Effect of the in-plane dissipation on the dynamics 
of Josephson vortices 

When a strong magnetic field is applied perpendicularly to 
the c-axis of BSCCO single crystal, the Josephson vortices are 
induced and form the triangular lattice. Driving by the Lorentz 



force induced by a transport current, the Josephson vortices 
move and excite Josephson plasma. The motion of Joseph- 
son vortex lattice provides an alternative routine to achieve 
a strong THz radiation. 32 Due to the strong inter-vorte x re- 
pulsion, the Josephson vortices favor the triangular lattice^E] 
and the radiation is weak. The rectangular lattice is observed 
in a small mesa 35 where the surface potential favors the rectan- 
gular lattice. Here we investigate the possible synchronization 
of the Josephson vortices by the in-plane dissipation. 

The simulation results are presented in Fig. [5] When the 
bias current increases, the Josephson vortices evolve toward a 



rectangular lattice. The rectangular lattice is achieved only by 
a high bias current, which is difficult to realize experimentally 
due to the strong self-heating effect. Our simulations are con- 
sistent with the analytical results obtained by Koshelev and 
Aranson 17 , who found that the rectangular lattice can only be 
stability in a high velocity region. The in-plane dissipation 
does not stabilize the in-phase plasma oscillation or rectangu- 
lar lattice of Josephson vortices due to the strong inter- vortex 
repulsion. The realization of the in-phase oscillation in the 
case of Josephson vortices is still an open problem and re- 
quires a new mechanism. 
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